Photoemission spectra of many-polaron systems 
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The cross over from low to high carrier densities in a many-polaron system is studied in the frame- 
work of the one-dimensional spinless Holstein model, using unbiased numerical methods. Combining 
a novel quantum Monte Carlo approach and exact diagonalization, accurate results for the single- 
particle spectrum and the electronic kinetic energy on fairly large systems are obtained. A detailed 
investigation of the quality of the Monte Carlo data is presented. In the physically most important 
adiabatic intermediate electron-phonon coupling regime, for which no analytical results are avail- 
able, we observe a dissociation of polarons with increasing band filling, leading to normal metallic 
behavior, while for parameters favoring small polarons, no such density-driven changes occur. The 
present work points towards the inadequacy of single-polaron theories for a number of polaronic 
materials such as the manganites. 
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I. INTRODUCTION 

In recent years, it has become widely accepted that 
electron-phonon (EP) interaction plays a crucial role in 
systems with strong electronic correlations, such as high- 
temperature superconducting cupratesi or colossal mag- 
netoresistive manganites.-- A large amount of the avail- 
able data on optical and transport properties has been 
interpreted in terms of polaronic states. In view of 
the complex underlying physics, however, the applica- 
tion of single-polaron theories seems to be problematic. 
On the other hand, a reliable theory of strongly cou- 
pled EP systems with finite carrier densities is not avail- 
able yet. While limiting cases, such as weak coupling 
or large phonon frequencies — compared to the electronic 
bandwidth — can to some degree be described theoreti- 
cally, especially the adiabatic regime of small phonon 
frequencies at intermediate-to-strong EP coupling con- 
stitutes a long-standing, open problem. Since analytical 
results can only be obtained at the expense of uncon- 
trolled approximations, progress toward an understand- 
ing of many-polaron systems requires unbiased numerical 
methods. Owing to the availability of high performance 
computers, the latter can yield insight into the properties 
of basic microscopic models such as the spinless Holstein 
model considered here. Nevertheless, apart from work on 
the half- filled band case, where the physics is dominated 
by a cross over from a Luttinger liquid to a Peierls state, 
the many-polaron problem has received little attention 
in the past. 

Motivated by this situation, in the present work, we 
address the important issue of the effect of carrier den- 
sity on the character of the quasiparticles of the system. 
While for very strong couplings no significant changes are 
expected due to the existence of rather independent small 
(self-trapped) polarons with negligible residual interac- 
tion, a density-driven cross over from a state with large 
polarons to a metal with weakly dressed electrons may 



occur in the intermediate coupling regime. This problem 
has recently been investigated experimentally in terms of 
optical measurements on La2y3(Sr/Ca)j^^3Mn03 filmsi^ 
While small polarons have been identified as the charge 
carriers in the low-temperature metallic phase of the Ca 
doped system, the weaker electron-lattice coupling in the 
Sr doped system gives rise to large-polaron-like behav- 
ior of the optical conductivity. The effect of the finite 
carrier density on such extended polaronic quasiparticles 
has been analyzed^ using a recently developed weak- 
coupling theoryi^ 

For the case of the Holstein model, the abovemen- 
tioned density-driven transition from small to large po- 
larons and finally to weakly EP-dressed electrons is ex- 
pected to be possible only in one dimension (ID), where 
large polarons extending over more than one lattice site 
exist at weak coupling. Existing work indicates that in 
higher dimensions, the electrons remain quasifree at weak 
coupling, and become self-trapped above a critical value 
of the EP coupling strengthiSiLSiS By contrast, approxi- 
mate variational calculations suggest that large polarons 
may also exist in D > l,iSiiiii^ The situation is different 
for Frohlich-type modelsi^ii^ with long-range EP inter- 
action, in which large polaron states exist even for strong 
coupling and in higher dimensions. 

Here, we make use of quantum Monte Carlo (QMC) 
and exact diagonalization (ED) methods, and calculate 
the photoemission spectra for finite clusters at zero and 
finite temperature. Due to the limitations of existing al- 
gorithms, we develop a new grand-canonical QMC ap- 
proach which is free of autocorrelations and therefore 
enables us to study small phonon frequencies and low 
temperatures. 

The organization of this paper is as follows. In Sec. ITU 
we present the model and discuss previous work, while 
Sec, mil contains details about the methods we employed. 
Section llVI gives a discussion of our numerical data, and 
finally, we conclude in Sec. The new QMC technique 
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is described in the Appendix. 



II. MODEL 

The spinless Holstein model is defined by the Hamil- 
tonian 

where c\ (c^) creates (annihilates) a spinless fermion at 
lattice site i, fcj (6^) creates (annihilates) a phonon of 
energy {h= 1) at site i, and n- = c|cj. While the first 
term describes hopping processes between neighboring 
lattice sites the second term corresponds to the 

elastic and kinetic energy of the phonons. Finally, the 
last term in Eq. ^ represents a local coupling of the 
lattice displacement to the electron occupation = 0, 1. 

The parameters of the model are the hopping integral 
t, the phonon energy Wq, and the coupling constant g. 
The latter is related to the polaron binding energy 
via g = Ep/ujQ. By defining the dimensionless quan- 
tities A — Ep/2tD, and ZUq = Wg/^' ^^^^ with two 
independent parameters. In the sequel, we express all 
energies in units of t. The lattice constant is taken to be 
unity, and periodic boundary conditions are applied. 

The Holstein Hamiltonian ^ represents a generic 
model for a many-polaron system. Due to the restric- 
tion to spinless fermions, on-site bipolaron formation is 
suppressed. 

Quantum fluctuations of the lattice enable the elec- 
trons to hop even in the strong-coupling regime, while 
in models with classical phonons translation symme- 
try is broken, so that no itinerant polaron states exist. 
Along this line, the half-filled spinless Holstein model has 
been studied in one dimension using QMCjiSiiLiS EDjiS 
DMRG^o and variational methods-SLSSiSiSiSi Similar to 
the spinfuU model at half filling, the spinless model with 
one electron per two lattice sites displays a quantum 
phase transition from a metallic state to a Peierls in- 
sulating state with charge-density-wave (CDW) order at 
T ~ 0. While the insulating CDW state of the spin- 
less model is stable for any A > in the adiabatic limit 
ZUq = 0, it can be destroyed by quantum phonon fiuctu- 
ations for tUg > 0^ Capone et al.^^ used ED to study 
the cross over from a single polaron to a many-electron 
system away from half filling. Although their work is ac- 
tually for the Holstein-Hubbard model, results are very 
similar to the spinless Holstein model in the limit of large 
U considered in Ref. 26. The authors found that the crit- 
ical coupling for polaron formation is unaffected by the 
electron density in the antiadiabatic regime. Finally, we 
would like to point out previous analytical efforts to take 
into account squeezing and correlation effects in many- 
polaron systems iSLSSiSi 



III. METHODS 

In this work we use QMC, and ED in combination with 
the kernel polynomial method (KPM), to calculate, in 
particular, single-particle spectra. These two approaches 
are characterized by different advantages and shortcom- 
ings. ED in the form used here yields exact spectra at 
T — with extremely high energy resolution. However, 
even for relatively small clusters of ten sites — the size 
used here — this method requires parallel supercomputers 
to handle the enormous Hilbert space (matrix dimension 
« 10^-10^°). Consequently, the momentum resolution, 
determined by the cluster size, is rather crude. 

In contrast, the grand-canonical QMC simulations can 
be performed on personal computers, and for larger clus- 
ters. The drawback is that (a) the method is limited to 
finite temperatures, and (b) the calculation of spectral 
properties requires the analytic continuation from imag- 
inary time to real frequencies, which is an ill-conditioned 
problem. Therefore, regularization techniques such as 
the maximum entropy method (MEM) used here have 
to be applied, which offer only limited energy resolution 
and may introduce uncontrolled errors. For this reason, 
we will also present a direct comparison of ED and QMC 
results. 



A. Quantum Monte Carlo 

In applying existing QMC methods for the Holstein 
model to the parameter regime considered here, one faces 
extremely long autocorrelations times, as well as large 
statistical errorsi^^ This makes it very difficult to ob- 
tain accurate results for spectral properties, due to the 
use of the MEM. To overcome these problems, the algo- 
rithm presented in the Appendix is based on the canon- 
ical Lang-Firsov transformatiouiS The latter explicitly 
contains the shift of the lattice equilibrium position in 
the presence of an electron, which is favorable for the 
QMC simulations. Moreover, the absence of a direct 
coupling between the electron density and the lattice dis- 
placements in the transformed Hamiltonian ljA2|l permits 
exact, uncorrelated sampling of the phonon degrees of 
freedom (momenta) in terms of principal components. 

Here we have mainly used a value oi fit = 8 for the 
inverse temperature /3 = 1/k^T. For the critical param- 
eters studied below, accurate simulations at even lower 
temperatures become very demanding due to a minus- 
sign problem (Sec. lA Sjl . Nevertheless, we shall present 
results for the density of states at (3t — 10 and half fill- 
ing but for a reduced cluster size. Note that already 
for /3t = 8 we have k^T/t = 0.125 < ZUg for the values 
(Ug > 0.4 used in the sequel, so that direct thermal exci- 
tation of phonons is expected to be rather unimportant. 

For the Trotter decomposition, we chose Ar = 0.05 
(see Appendix) for static observables. Owing to the 
larger numerical effort associated with the calculation of 
dynamical properties, for the latter we use At = 0.1 — 
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the larger error being small compared to the uncertainties 
introduced by the MEM inversion. The Trotter error in 
the Holstein model becomes important mainly for large 
phonon frequencies and very strong EP interaction. 
Since we do not consider this regime here, the abovcmcn- 
tioned values of At yield satisfactorily small systematic 
errors. 

Away from half filling, the electron density has to be 
adjusted to the values of interest by varying the chemical 
potential Due to the computational effort associated 
with this trial-and-error procedure, the actual band fill- 
ing n', say, is usually only very close to a desired value, 
e.g., n = 0.4, but not exactly the same. Furthermore, 
within the QMC simulations, the filling can only be de- 
termined up to a statistical error An', similar to other 
observables. For the results presented here, the relative 
deviation {n — n') /n of the actual value in the QMC sim- 
ulations from the value n reported is always smaller than 
5 X 10^'^. The same is true for the relative statistical 
error An' /n. 

To perform the continuation of the one-electron Green 
function from imaginary time to real frequencies, we 
use the MEM (for a review see, e.g., Ref. El). The 
algorithm employed here incorporates an improved ev- 
idence approximation,'^^ which is less susceptible to 
overfitting — resulting in peaks whose position, weights 
and even existence are not clear — than the conventional 
evidence approximation^ On the other hand, groups of 
narrowly spaced peaks are represented by an envelope if 
the accuracy of the QMC data is not sufficient. 

The proposal of an alternative evidence 
approximation — leading in general to stronger reg- 
ularization and, consequently, lower energy resolution 
than the conventional approximation — may not seem 
advantageous at first. However, considering the great 
variety of different spectra compatible with the QMC 
data within the error bars due to the extremely ill- 
conditioned analytical continuation problem, we argue 
that it rather corresponds to choosing the "most un- 
informative" solution in the spirit of the maximum 
entropy principle. In fact, we found that the alternative 
evidence approximation yields much more consistent 
results than the "classic" evidence approximation;^ 
when we subjected artificial test data to different noise 
while keeping the signal-to-noise ratio fixed. A detailed 
comparison of different regularization schemes available 
will be presented elsewhere. 

In order to scrutinize the resulting spectra, the follow- 
ing test has been performed. Following the idea of Silver 
et al.,^^ the results from the MEM applied to the orig- 
inal QMC data G{k, r) have been taken as an "exact" 
spectral function A{k, u), from which mock data G{k, t) 
for the imaginary-time Green function has been created 
by using Eq. HA19|) and adding Gaussian random noise 
of absolute size 10~'^-10~*, comparable to the statisti- 
cal errors of the QMC results. Again using the MEM 
to invert these data, the thereby obtained spectral func- 
tion A{k,uj) allows to check the stability of the MEM 



result A{k,u!) with respect to statistical errors. We find 
that the features in our results for A{k, uj) remain sta- 
ble if G{k,uj) varies within the statistical errors, so that 
the data presented below are satisfactorily reliable apart 
from the aforementioned broadening and the limited en- 
ergy resolution. 

While our method is completely free of autocorrela- 
tions between subsequent phonon configurations, results 
for the imaginary-time Green function at different times 
are usually not statistically independent.'^^ However, in 
the present case, owing to the exact sampling and the 
large number of measurements, these correlations are 
very small and have been neglected since they have neg- 
ligible influence on the results. An additional, stringent 
test of the QMC results consists of a comparison to ED 
data which will be presented in Sec. IIVBI 

Finally, we note that all QMC simulations have been 
run on personal computers (e.g., Intel XEON 2600 MHz), 
with CPU times varying between several hours and a 
couple of days. 



B. Exact diagonalization 

The combination of ED and KPM has been described 
in detail in Ref.lSSl In this work, we are mainly interested 
in the single-particle spectral functions 

A+(fc,c.) = ;^i(<r^)|ci._j<))P 

x^[.-(i^M 
A-{k,.) = Y^\{¥,f-'^\c^^,\¥Z^^)\^ 

(2) 

where |^'|^°'') denotes the Ith eigenstate with iV, elec- 
trons, momentum k and energy e11°\ The cluster size 
has been fixed at N — 10, with as many as 15 dy- 
namical phonons and resulting truncation errors < 10~^ 
for the ground-state energy. "^^ Additionally, the conver- 
gence with the number of phonons has been monitored 
by means of the phonon distribution functioni^ It is im- 
portant to point out that such a small number of dy- 
namical phonons only yields accurate results for the pa- 
rameters considered here if the symmetric phonon mode 
with q — Ois separated from the Hamiltonian analytically 
before the diagonalization. '^^''^^ In the present case, this 
phonon mode — corresponding to a homogeneous shift of 
the oscillator equilibrium positions — together with that 
for q = TT represents the dominant contribution to phonon 
excitations, so that the separation is very efficient at re- 
ducing the computational effort. The largest matrix di- 
mension in the present calculations was about 3 x 10*, 
running on six compute nodes of a Hitachi SR8000-F1. 
The energy resolution of the KPM spectral function with 
512 moments has been enhanced by a factor of 8 by the 
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MEM)2Si resulting in a significantly better resolution com- 
pared to, e.g., the spectral decoding method. We would 
like to point out that in contrast to the case of QMC, the 
ED data to which the MEM is applied does not have a 
statistical error. 



IV. RESULTS 

It is long known that within the Holstein model Q), 
a single electron undergoes a cross over to a small po- 
laron at a critical value of the EP coupling>& The im- 
portant question to be addressed here is what happens 
if the density of electrons (or polarons) is increased. We 
do not expect significant changes in the strong-coupling 
regime due to the rather "localized" nature of the small 
polaron state. On the contrary, for intermediate EP in- 
teraction and ZZJq < 1, each electron is surrounded by a 
phonon cloud extending over several lattice sites. Now, 
at larger densities, a substantial overlap of the single- 
particle wavefunctions occurs, leading to a dissociation 
of the individual polarons and finally to a restructur- 
ing of the whole many-particle ground state. Further- 
more, it has been shown that squeezing effects reduce the 
polaronic band-narrowing in the many-electron case^Si 
thereby favoring such a scenario. Since within the Hol- 
stein model, a large polaron state is expected to be stable 
in ID, we will only consider this case in the sequel. 

In view of the situation in most polaronic materials, we 
mainly consider the adiabatic regime, taking oj^ = 0.4. 
To illustrate the important differences between oJg <C 1 
and 3> 1, however, some results for the antiadiabatic 
value ZUq = 4 will also be shown. 



A. Limiting cases 

As we shall see below, the spectra in the intermediate 
coupling region turn out to have a fairly complex struc- 
ture. It is therefore helpful to begin with the limiting 
cases of weak and strong EP coupling. 



1. Weak coupling 

For weak coupling A = 0.1, the sign problem is not 
severe (Sec. lA 3|l and the QMC simulations can easily 
be performed for large lattices with = 32, thereby 
making the dispersion of quasiparticle features well visi- 
ble. As indicated in Eq. (|A19I) . the MEM inversion yields 
A{k,ui — jj), where /i denotes the chemical potential. 

Figure n shows the evolution of the one-electron spec- 
tral function A{k^LO — /i) with increasing electron den- 
sity n. At first sight, we see that the spectra bear a 
close resemblance to the free-electron case, i.e., there is 
a strongly dispersive band running from —2t to 2t. The 
latter can be attributed to weakly dressed electrons with 
an effective mass similar to the noninteracting case. As 
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FIG. 1: One-electron spectral function A{k,uj~fj) from QMC 
for different band fillings n, N = 32, pt = 8, ujq — 0.4, and 
weak coupling A — 0.1. 
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FIG. 2: One-electron density of states p{uj ~ fi) from QMC 
for different band fillings n, N — 32, fit = 8, ZJo — 0.4 and 
A = 0.1. 



expected, the height (width) of the peaks increases (de- 
creases) significantly in the vicinity of the Fermi momen- 
tum fcp, which is determined by the crossing of the band 
with the chemical potential. 

Since the peaks in Fig. ^ are generally very sharp, it is 
important to check the results from the inversion with the 
MEM. By computing the approximate number of elec- 
trons corresponding to the fillings n = 0.1-0.4 of the 
32-site system, we can roughly determine k^. Although 
the nonzero EP interaction will affect the band dispersion 
E{k), the weak coupling considered here justifies the use 
of a rigid band approximation, as reflected in the free- 
particle character of the spectra in Fig. ^ 

For n = 0.1, i.e., about three electrons for N = 32, the 
Fermi momentum is expected to lie between k = tt/IQ 
and k = n/8, corresponding to the second and third curve 
from the top in Fig.^a). In fact, we see that the largest 
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FIG. 3: One-electron spectral function w — /x) from QMC 
for different band fillings n, N = 32, l3t = 8, ujq = 0.4, and 
strong coupling A = 2.0. 



peaks are those for k = n/lQ and tt/8, having almost 
the same height. Above and below these momenta, the 
peaks are noticeably smaller. A similar analysis carried 
out for the other fillings in Figs, ^b) -(d) also shows a 
very good agreement between the location of the highest 
peaks and the estimated values of fcp. 

Finally, we would like to point out that the apparent 
absence of any phonon signatures in Fig. ^ is not a de- 
fect of the MEM, but results form the large scale of the 
z-axis chosen. As a consequence, the peaks running close 
to the bare band dominate the spectra and suppress any 
small phonon peaks present. If we enlarge the resolu- 
tion, for all densities n = 0.1-0.4, we observe the band 
flattening2Si^2ii£ at large wavevectors which originates 
from the intersection of the approximately free-electron 
dispersion with the optical phonon energy occurring at 
w — ^ = Wq. 

To complete our discussion of the weak-coupling 
regime, we show in Fig. |21 the one-electron density of 
states (DOS) p{uj — fi) given by Eq. HA22|I . Clearly, for 
small n, there is a peak with large spectral weight at the 
Fermi level. In contrast, for large n, the tendency to- 
ward formation of a Peierls- (band-) insulating state at 
n = 0.5 suppresses the DOS at the Fermi level, although 
we are well below the critical value of A at which the cross 
over to the insulating state takes place at T = Oi^ The 
additional small features separated from fi by the bare 
phonon energy cJq will be discussed in the intermediate 
coupling case below. 



2. Strong coupling 

We now turn our attention to the opposite limit of the 
spinless Holstein model, namely the case of strong EP 
coupling A = 2. While calculations with many other nu- 
merical methods such as, e.g., standard grand-canonical 
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FIG. 4: Dispersion of the polaron band in A{k,uj — ji) for 
n = 0.05, A'' = 20 and /3i = 10 (drawn to a larger scale as 
compared to Fig. 3). All other parameters as in Fig. |21 



QMC or ED, become very difficult in this regime, the 
QMC approach used here yields quite accurate results 
even for large systems. The improved statistics in the 
strong-coupling regime is also reflected in the dependence 
of the average sign on A (cf. the discussion of Fig. IT^ 
in the Appendix). As before, we study the dependence 
of the spectral function on the band filling, again for 
N — 32. The results are depicted in Fig.|3| 

In the low-density regime n = 0.1 [Fig.|3Ja)], we expect 
the well-known almost flat polaron band having expo- 
nentially reduced spectral weight (given by in the 
single-electron, strong-coupling limit) which, neverthe- 
less, can give rise to coherent transport at T = 0. Unfor- 
tunately, the QMC/MEM has difficulties resolving such 
extremely weak signatures, especially at larger n, where 
the statistical errors are more noticeable, and for temper- 
atures such that the thermal energy is comparable to the 
polaron band width. Therefore the polaron band can 
hardly be resolved in Figs.Ota) or (b), and is not even 
detectable in Figs.EIc) and (d). However, if we consider 
a density n — 0.05, corresponding to a single electron on 
a cluster with = 20, at a lower temperature fit = 10, 
the QMC results do indeed give a polaron band around 
jjL (Fig. 21). Due to the sign problem, similar calculations 
for n > 0.1 would be very demanding. 

Besides, the spectrum consists of two incoherent fea- 
tures located above and below the chemical potential, 
which reflect the phonon-mediated transitions to high- 
energy electron states. Here, the maximum of the pho- 
toemission spectra (w — /i > 0) follows a tight-binding 
cosine dispersion. The incoherent part of the spectra is 
broadened according to the phonon distributions. We 
would like to mention that the spectrum of a single po- 
laron at strong EP coupling has been studied analytically 
by Alexandrov and Ranninger.'*^ 

With increasing band filling [see panels (b) and (c)], 
the chemical potential is shifted to lower energies cx Epn, 
but is expected to be still located in a narrow band with 
little spectral weight. Regardless of the aforementioned 
shortcomings, the QMC results provide a clear picture of 
the physics. There exists a finite gap to the photoemis- 
sion (inverse photoemission) parts of the spectrum, so 
that the system typifies as a polaronic metal. We shall 
see below that a completely different behavior is observed 
at intermediate coupling. Notice that the incoherent in- 
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FIG. 5: One-electron density of states p{u — ^) from QMC for 
different band fillings n and cluster sizes A'', (31 = — 0.4 

and A = 2.0. 

verse photoemission (photoemission) signatures are more 
pronounced at small (large) wavevectors. 

Finally, for n — 0.4 [Fig. I^Jd)], the incoherent features 
lie rather close to the Fermi level, thus being accessi- 
ble by low-energy excitations. Now, the photoemission 
spectrum for fc < 7r/2 becomes almost symmetric to the 
inverse photoemission spectrum for fc > 7r/2 and already 
reveals the gapped structure which occurs at n = 0.5 due 
to the CDW-formation accompanied by a Peierls distor- 
tion. 

As in the weak-coupling case discussed above, the 
properties of the system also manifest itself in the DOS, 
shown in Fig.|Sl Owing to the strong EP interaction, the 
spectral weight at the chemical potential is exponentially 
small for all fillings n. 

At half filling, the DOS exhibits particle-hole symme- 
try, and the system can be described as a Peierls insula- 
tor, consisting of a polaronic superlattice. In contrast to 
the weak-coupling case, the ground state is characterized 
as a polaronic insulator rather than a band insulator. 
Moreover, here the gap at the Fermi level persists even 
at the finite temperature used in the simulation, while 
for A = 0.1 and n = 0.5 the insulating state found at 
T = is destroyed by thermal fiuctuations (Fig. [21 . Note 
that, owing to the sign problem, the result for n — 0.5 
in Fig. [S] is for a reduced cluster size N — 20. However, 
in the strong-coupling regime considered here, finite-size 
effects are expected to be very small. 

B. Intermediate coupling 

As discussed in the introduction, a cross over from a 
polaronic state to a system with weakly dressed electrons 
can be expected in the intermediate coupling regime. 
Here we choose A = 1, which corresponds to the critical 
value for the small-polaron cross over in the one-electron 
problem (cf., e.g.. Fig. 5 in Ref. l2l). 



1. Quantum Monte Carlo 

We first discuss the QMC results. Owing to the sign 
problem, which is particularly noticeable for A = 1 (see 
Appendix), we have to decrease the system size as we 
increase the electron density n. 

Figure |H1 shows the spectral function for A = 1 and 
increasing band filling. Owing to the overlap of large po- 
larons in the intermediate coupling regime, we start with 
the very low density case n — 0.05 [Fig.|B|[a)]. Compared 
to the behavior for A = 2 [Fig.|3{a)], we notice that the 
polaron band now lies much closer to the incoherent fea- 
tures, and that there is a mixing of these two parts of the 
spectrum at small values of fc. Nevertheless, the almost 
flat polaron band is well visible for large k. 

With increasing density, the polaron band merges with 
the incoherent peaks at higher energies. This is the 
above-anticipated density-driven cross over from a pola- 
ronic to a (diffusive) metallic state, with the main band 
crossing the Fermi level. At n = 0.4, we find a band with 
large effective width ranging from oj « —At to w « 2t. A 
more detailed discussion of the spectral function will be 
given below when we present ED results. 

Further information about the density dependence can 
be obtained from the one-electron DOS. The latter is 
presented in Fig. [7| for different fillings n — 0.05-0.5. 
As in Fig. |S1 the cluster size is reduced with increasing 
n in order to cope with the sign problem. To illustrate 
the rather small influence of finite-size effects, Fig.[71also 
contains results for N — 10. 

For low density n — 0.05, the DOS in Fig.[7|lies in be- 
tween the results for weak and strong coupling discussed 
above. Although the spectral weight at the chemical po- 
tential is strongly reduced compared to A = 0.1, p{0) is 
still significantly larger than for A — 2.0. 

When the density is increased to n = 0.2, the DOS 
at the chemical potential increases, as a result of the 
dissociation of the polarons. As we increase n further, a 
pseudogap begins to form at /i, which is a precursor of 
the CDW gap at half filling and zero temperature. 

Finally, in the case of half filling n = 0.5, the DOS 
has become symmetric with respect to fj,. There are 
broad features located either side of the chemical po- 
tential, which take on maxima close to w — /i = ±Ep. 
However, apart from the strong-coupling case, where the 
single-polaron binding energy is still a relevant energy 
scale, the position of these peaks is rather determined by 
the energy of the upper and lower bands, split by the for- 
mation of a Peierls state. The gap of size ~ A expected 
for the insulating charge-ordered state at T = is par- 
tially filled in due to the finite temperature considered 
here. 

Furthermore, we find additional, much smaller features 
roughly separated from /x by the bare phonon frequency 
(jjg, whose height decreases with decreasing temperature, 
as revealed by the results for f3t = 10 (Fig. [7J|. These 
peaks — not present at T = (see Ref. l37j) — arise from 
thermally activated transitions to states with additional 
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FIG. 6: One-electron spectral function A{k, u) - 
and intermediate coupling A = 1.0. 



' fj,) from QMC for different band fillings n and cluster sizes A*', /3t = 8, ujo ~ 0.4, 
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FIG. 7: (color online) One-electron density of states piuj — fi) 
from QMC for different band fillings n, cluster sizes TV and 
inverse temperatures /3. Here tJo = 0.4 and A — 1.0. 



phonons excited, and are also visible in Figs. |2] and |5l 
While for weak coupling [Ep = 0.2, A = 0.1, Fig. El, the 
maximum of these features is almost exactly located at 



couphng [Ep = 2, A = 1, Fig. [7], and finally to \uj — 
/i| « 2.5wq for strong couphng [Ep — 4, X = 2, Fig. |5]. 
Although the exact positions of the peaks are subject to 
uncertainties due to the MEM, this reflects the shift of 
the maximum in the phonon distribution function with 
increasing coupling. The MEM yields an envelope of the 
multiple peaks separated by Wg. 

At this point, we would like to mention the exact 
relation^ 



M,{k) 



duj u) A{k, Lo) 



2Epn 



(3) 



for the first moment of the one-particle spectral func- 



tion of the spinless Holstein model, with e^. 



-2t cos k. 



-Ml 



, it moves to 



-/i| w 1.25ciJg for intermediate 



While the zeroth moment is identical to the normaliza- 
tion of A{k,uj) for each k, M^{k) depends in a nontriv- 
ial way on the parameters of the system. In principle, 
Eq. Q may be used as a boundary condition in the 
MEM.'*'^ However, in the present case, we have found 
that while the normalization is virtually exact for all re- 
sults shown, the first moment deviates from the exact 
values. Moreover, the use of the exact values for 
as a condition in the inversion causes the MEM not to 
converge for some parameters. Additional calculations 
on small systems reveal that this disagreement results 




FIG. 8: (color online) One-electron spectral functions A~{k,u!) (red dot-dashed lines) and A^{k,Ltj) (black thin solid lines) 
from ED for different band fillings n, N = 10, ZJo ~ 0.4, and intermediate coupling A — 1.0. Also shown are the integrated 
spectral weights of A~ and A'^ (dashed lines), as well as A{k,uj) from QMC for A*' = 10 and (31 = 9, (blue thick solid lines). 



from the finite Trotter error, despite the relatively small 
value At = 0.1 used. Away from parameters such that 
K, 0, as determined from the MEM spectra devi- 
ates from the exact values by maximally 10-20 percent. 
Nevertheless, the peak positions and weights in the spec- 
tra shown are fairly accurate, as illustrated by the com- 
parison with ED below. Finally, we would like to point 
out that calculations of A(fe, uj) for Ar < 0.1 become very 
time-consuming due to the increasing requirement in disk 
space (see Appendix) and longer simulation times. 

Despite the wealth of information contained in the 
QMC results presented so far, it is obvious that an inter- 
pretation of the various excitations is far from being easy. 
In addition to the finite temperature, the QMC method 
used here does not allow to separately calculate the pho- 
toemission and inverse photoemission parts A~ and , 
respectively. Moreover, the use of the MEM limits the 



energy resolution and introduces uncertainties concern- 
ing the positions and weights of structures in the spectra. 
These circumstances make it difficult, e.g., to distinguish 
between phononic and electronic contributions, especially 
if A « ZZJg. To gain further insight, wc therefore supple- 
ment the QMC data with ED results. 



2. Exact diagonalization 

In order to overcome the abovementioned limitations of 
QMC, here we use ED in combination with the KPM, as 
described in Sec. lIIIBl to calculate the one-particle spec- 
tral functions defined by Eq. (0). To test the quality of 
the QMC spectra, we directly compare the two methods, 
keeping in mind that the QMC calculations have been 
performed at a finite temperature, whereas ED yields 
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CO/ t 

FIG. 9: (color online) Comparison of the integrated spectral 
weight obtained from the QMC and ED data of Fig. |HIa) for 
three different momenta k. 



ground-state results. 

Figure El displays the ED results for the same parame- 
ters as in Fig. El i.e., cZJq = 0.4, A = 1, and n = 0.1-0.4. 
Additionally shown are the integrated spectral weights 
of A~ and yl+ (dashed lines), as well as QMC results for 
the same system size and pt — 8 (solid thick lines), which 
have been shifted by /i. The Fermi level lies between the 
highest-energy peak in and the lowest-energy peak in 
A+. 

We see from Fig.jHlthat the QMC results are compati- 
ble with the ED data. Because of the finite temperature 
and the use of the MEM, QMC, of course, cannot repro- 
duce the sharp ED peaks. To further test the agreement 
between the two methods, we compare in Fig. the in- 
tegrated spectral weight of the data shown in Fig. |HIa), 
and find a very satisfactory accord. Moreover, this agree- 
ment becomes even better for larger densities n, for which 
the sharp peaks in the ED data evolve into bands [see 
Fig. Ob) -(d)]. 

In addition to the possibility of distinguishing between 
A~^ and A~ — crucial for an identification of polaron 
bands at intermediate coupling — the higher energy reso- 
lution of ED also allows to resolve signatures of phonon 
excitations. 

Starting with the case n = 0.1 [Fig. jSJa)], we notice 
that there is a polaron band at the Fermi level, exci- 
tations to/from which are given by the highest (lowest) 
peak of A" (A+). The photoemission part A~ reflects 
the Poisson-like phonon distribution of the one-electron 
ground state {n — 0.1 and N — 10). The integrated spec- 
tral weight gives a measure for the electronic weight of 
the various poles in the one-electron spectrum. For ex- 
ample, it reveals that the phonon peaks in A^ have very 
little weight. In fact, the integrated weight jumps to a fi- 
nite value at the first peaks near the Fermi level for k = 
and k — 7r/5, while it changes very little as one moves 
further down in energy. For k > 7r/5, the small spec- 
tral weight contained in A~ is continuously distributed 
among the phonon peaks. 

We also observe the well-known flattening of the po- 
laron band at large values of k. Similar to the one- 
electron case^^ the low-energy states have mostly elec- 
tronic character at small fc, and become mostly phononic 
at large values of k. While this effect is expected to occur 
in the low-density regime, we find that is persists even 
for n = 0.4 [Fig.||ld)]. Moreover, Fig. reveals that the 
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FIG. 10: (color online) One-electron spectral functions 
A~{k,Lu) and A+{k,uj) from ED for TV = 10, ZJo = 4.0, and 
A — 2.0. Here n = 0.1 {A~: red dot-dashed lines, A^: black 
solid lines) and n = 0.3 {A~: green dashed two-dotted lines, 
A^: blue dotted lines). Dashed lines correspond to the inte- 
grated spectral weight of A~{k,uj) and A'^{k,uj), respectively. 



maximum in the incoherent contribution follows closely 
the free-electron dispersion for all densities n — 0.1-0.4. 

With increasing electron density, the equally spaced 
peaks in A" broaden significantly, until at n = 0.4 they 
have merged to form a broad band. The polaron band 
is still visible at n = 0.2 [Fig. |HJb)], but becomes in- 
distinguishable from the incoherent excitations at even 
larger n. Eventually, at n = 0.4 [Fig. Od)], we find 
an almost symmetric behavior of A~^ and A~ with re- 
spect to the Fermi level at k — 27r/5. Compared to the 
strong-coupling case [Fig. Old)], there exist incoherent 
excitations with almost zero energy. This indicates the 
metallic character of the system expected for large po- 
laron densities and intermediate EP coupling. 



C. Antiadiabatic regime 

Until here, we have only presented results for the adi- 
abatic case ZUq — 0.4. Although the latter is most rele- 
vant for strongly correlated materials, it is important to 
compare the above findings to the antiadiabatic strong- 
coupling regime. To this end, we chose ZUp = 4 as well 
as A = 2. The condition for the formation of a single 
small polaron in this case, however, is g > 1, because the 
phonon frequency is the important energy scale. In con- 
trast, the criterion A > 1 is purely based on the balance 
of kinetic and potential energyi^24^ii^ 

The spectral function for n = 0.1 and n — 0.3 is dis- 
played in Fig. 1101 for 5=1, i.e., at the transition point 
where self-trapping occurs for a single particle. Contrary 
to the case A = 1 and oJq = 0.4 considered above, the 
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FIG. 11: Normalized kinetic energy as a function of EP cou- 
pling A for different band fillings n and adiabatic ratios ZZJo- 
Also shown are the results for a single polaron, obtained with 
the method of Ref.l23. Errorbars are smaller than the symbols 
used and Ar = 0.05. 

absorption features are now separated by the phonon en- 
ergy Wq, with their width being determined by electronic 
excitations. 

A comparison of the spectral functions for n = 0.1 and 
n = 0.3 in Fig. ^1 indicates that there is no density- 
driven cross over of the system as observed in the adi- 
abatic case. In particular, owing to the large phonon 
energy, there are no low-energy excitations close to the 
polaron band, so that the latter remains well separated 
from the incoherent features even for n = 0.3. Further- 
more, the spectral weight of the polaron band also re- 
mains almost unchanged as we increase the density from 
n = 0.1 to n = 0.3. Consequently, almost independent 
small polarons are formed also at finite electron densities, 
in accordance with previous findings of Capone et alm^ 
for small systems. 

D. Kinetic energy 

Previously, the electronic kinetic energy has been stud- 
ied extensively to monitor the polaron and bipolaron 
cross over in the Holstein model with one and two elec- 
trons (see, e.g., Refs-lolandli^. Here, at least in the adia- 
batic regime, polaron formation manifests itself in a drop 
of the kinetic energy around the critical coupling. In the 
following, however, we illustrate that the kinetic energy 
[Eq. (|A16|l ] does not refiect the significant changes in the 
one-particle spectrum we found with increasing density 
in the intermediate couphng regime (see Sec. lIVBl) . 

To compare different band fillings, we consider the nor- 
malized quantity E^(n, X)/ E^^{n,0), shown in Fig. II II as 
a function of EP coupling for different fillings and phonon 
frequencies, for a fixed system size = 16 and pt = 8. 

In the adiabatic regime cUq = 0.4, owing to the overlap 
of the phonon clouds surrounding the individual quasi- 
particles, the carriers clearly become more mobile with 



increasing density, but the change in the kinetic energy is 
nevertheless surprisingly small. In addition to these two 
fillings, we have also included in Fig. ^2 the result for a 
single polaron with the same parameters, calculated with 
the QMC method of Rcf. which is indistinguishable 
from the result for n = 0.1. 

As expected from the almost identical spectra for n = 
0.1 and 0.3 presented in Sec. IIV CI for the antiadiabatic 
regime, is almost independent of n in this case. 



V. CONCLUSIONS 

Strongly correlated electron-phonon systems in depen- 
dence on the carrier density have been studied by means 
of unbiased numerical methods. Using exact diagonaliza- 
tion together with the kernel polynomial method, as well 
as a novel quantum Monte Carlo approach, we were able 
to present accurate results for the one-dimensional spin- 
less Holstein model in the most difficult regime of small 
phonon frequencies and intermediate to strong electron- 
phonon interaction. 

We have calculated the one-particle spectral function, 
the density of states and the electronic kinetic energy 
at different band fillings. While exact diagonalization 
is restricted to relatively small clusters, thereby limit- 
ing the momentum resolution of the spectra, quantum 
Monte Carlo can be used to study larger systems, but 
yields a lower energy resolution for dynamical quantities. 
The reliability of the results from the maximum entropy 
inversion has been scrutinized, and we find a very good 
agreement between the two methods. 

In the adiabatic case, characteristic for, e.g., the man- 
ganites, we observe a carrier density-driven cross over 
from a polaronic state to a metallic system at weak and 
intermediate electron-phonon coupling. Note that the 
nature of the quasiparticles is changed in the course of 
this transition. This points toward an insufficiency of 
simple single-polaron theories to explain experimental re- 
sults for polaronic materials. In fact, very similar effects 
due to carrier interaction have recently been observed 
experimentally for the manganites.'^'^ 

On the contrary, for large phonon frequencies and 
strong coupling, the individual polarons remain virtually 
unaffected upon increasing the electron density. Conse- 
quently, the system remains polaronic even at large band 
fillings. Finally, the kinetic energy is shown to be too 
insensitive in order to reveal the above cross over. 

For weak electron-phonon coupling, the system be- 
comes a Peierls or band insulator at half filling, whereas 
for strong coupling a polaronic insulating state with an 
exponentially small polaron band at the Fermi level, and 
a finite gap to the incoherent excitations, prevails for all 
fillings. 

Finally, reliable calculations of the optical conductivity 
are highly desirable in order to explain the large amount 
of available experimental data for polaronic materials. 
Work along this line is in progress. 
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APPENDIX A: QUANTUM MONTE CARLO 

In this Appendix we extend the one-electron QMC al- 
gorithm developed in Ref. [2^ henceforth also referred to 
as I, to the spinless Holstein model with many electrons. 
We begin by applying the Lang-Firsov transformation^^ 
v = exp(i7 n^p^) with 7^ = 4AD/iUQ to the Hamil- 
tonian ^ in first quantization and with dimensionless 
phonon variables x and p, which takes the form^^ 



H 



-^E^k + 2 



(Al) 



The coupling constant a is related to that in Eq. by 
a — \/2gLUg. Using the result H = vHv^ (see I), we 
define the grand-canonical Hamiltonian 



n 



n(Pi-Pj) 



K 



P=P„+P^ 



(A2) 



where fi denotes the chemical potential and Ep is the 
polaron binding energy. 

For the case of a half-filled band, the chemical poten- 
tial is given by /i = —Ep. Away from n = 0.5, fi has 
to be adjusted accordingly to yield the right density of 
electrons n. 

We would like to point out that in the spinfuU case, 
the transformed Hamiltonian would contain an attractive 
Hubbard term. For the determinant QMC method to be 
applicable, the latter has to be decoupled using auxiliary 
fields,'*^ In contrast, for the spinless model considered 
here, the algorithm is almost identical to the one-electron 
case. The following derivation assumes a hypercubic lat- 
tice in D dimensions, consisting of sites. 



We use the Suzuki- Trotter decomposition^ 

L 



(A3) 



with At = P/L. The trace appearing in the partition 
function Z = Tr e~^'^ can be split up into a bosonic and 
a fermionic component leading to 



xn(p. 



AtK -AtP„-AtP^-AtI 



i=l 



where dp^ = dp^ ^, and with periodic boundary con- 
ditions Pj^^i = Pi- We have chosen the time ordering, 
which is arbitrary, in increasing order. The phonon co- 
ordinates X in Zj^ can be integrated out analytically in 
the same manner as in I. Moreover, the momenta p can 
be replaced by their eigenvalues on each time slice, and 
the partition function takes the form 



(A4) 



2l^CJ Vpe-'^^'-Tv.iB.B^ ---B^) 
with Vp = dp^, C = [27r/(tJoAT)]^-^/2 and 



=i7(Pi,x-Pi.x) 



The bosonic action S-^ may be expressed in terms of prin- 
cipal components (see I) . The fermion degrees of freedom 
can be integrated out exactly leading tc^* 



Trf(B,...i?J = det(l + i?, ••• i?J 
= dct(l + f7), 

where the x matrix B^ is given by 
B^ = D^kDIV 

T T T 

with 



(A5) 
(A6) 



Here h^^ is the usual tight-binding hopping matrixes 
To save computer time, we employ the checkerboard 
breakup 



(A7) 



Using Eq. (|A7I) . the numerical effort scales as N'^^ in- 
stead of N^^ , while the error due to this additional ap- 
proximation is of the same order (Ar)^ as the Trotter 
error in Eq. ||A3|I . 
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Defining the bosonic and fermionic weights Wy^ 



-ArSb 



and — det(l + f2) respectively, the partition 



function finally becomes 



C / Vp t«. Wf 



(A8) 



One of the advantages of Blankenbecler et aZ.'s*^^ formal- 
ism as well as the current approach is the close relation 
to the one-electron Green function 



Working in real space and imaginary time, we hav o^^i^^ 
G1,^{crc])^{l + n)T^\ (AlO) 

and 



G".. = 5-. 



(All) 



We would like to mention that despite the formal sim- 
ilarity to the method of Blankenbecler et al.^ the nu- 
merical realization of the present approach is quite dif- 
ferent. While the grand-canonical method of Ref. 
benefits enormously with respect to performance from 
a local updating scheme for the phonons, here we use 
a global updating in terms of the principal components 
together with the reweighting method described in I. Al- 
though this requires us to recalculate the full matrix 
in each MC step, the resulting statistically independent 
configurations clearly outweigh the loss in performance, 
especially for small uj^ for which autocorrelation times 
can exceed 10^ sweeps when using the original grand- 
canonical QMC methodi^ The approach proposed here 
allows one to perform uncorrelated measurements after 
each update. Moreover, our calculations show that nu- 
merical stabilization by means of a time-consuming sin- 
gular value decomposition is not necessary for all param- 
eters considered in Sec. IIVI Finally, owing to the phase 
factors in the transformed hopping term, all the matrices 
become complex-valued. An important practical test in 
simulations therefore is the reality of expectation values 
of observables, i.e., the disappearance of the imaginary 
part within statistical errors. 

The phase factors in the hopping term [Eq. ljA2p ] give 
rise to a minus-sign problem, i.e., Wf can take on negative 
values, whereas Wy^ is always positive. While for a fixed 
number of electrons the sign problem diminishes with 
increasing system size,^^ the opposite is true if we fix the 
electron density. In that case, we find an exponential 
decrease of the average sign with increasing cluster size, 
and a strong dependence on the band filling fSec. lA 3ll . 



2. Observables 

The calculation of observables within the formalism 
presented here is similar to the standard determinant 



QMC method An important difference is that we 
have to use the Lang-Firsov-transformed operators, i.e.. 



(O) --Tr(Oe-'^^). 



(A12) 



As mentioned above, the MC sampling is based purely on 
the bosonic weight w^^. This corresponds to a reweight- 
ing of the probability distribution so that the fermionic 
weight Wf is treated as part of the observables As for a 
single electron,^- the overlap of the probability distribu- 
tions involved is large enough to permit this procedure. 
Hence, we have 



(0) = 



'f/b 



{w. 



(A13) 



f/b 



where the expectation value {0)y^ for an equal-time ob- 
servable O is defined as 



a. Static quantities 



(A14) 



We begin with the electron density 



(A15) 



The expectation value in Eq. I|A15|I can be calculated 
from the diagonal elements of the Green function 
[Eq. (|ITO|) ]. i.e., {fi^) = {G\^). Similarly, the absolute 
value of the electronic kinetic energy per site is given by 



(A16) 



Equal-time two-particle correlation functions such as 
Pi^) — 'Ylii{''^i^i+s) n^^-y be calculated with the help of 
Wick's TheoremiSia 



h. Dynamic quantities 

An observable of great interest is the time-dependent 
one-particle Green function 



which is related to the spectral function 

1_ 



through 

G\k,T)-. 



duj ■ 



-lniG\k,u} - ^i) 



1 + e-/3(i^-/j) 



(A17) 



(A18) 



(A19) 
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tigation of the dependence of the average sign on the 
various parameters. 

While Hamiltonian l)A2|l is symmetric with respect to 
a particle- hole transformation for /z = — -Ep (half filling), 
this symmetry is broken if we use the checkerboard ap- 
proximation l|A7|l . so that n ^ 0.5 for the above choice 
of the chemical potential. To simplify calculations, the 
results for the average sign at n = 0.5 presented below 
have therefore been obtained using the exact hopping 
term. For general band fillings, we find that the checker- 
board breakup requires different values of /x to perform 
simulations at the same electron densities, but results are 
identical within statistical errors if /i is adjusted accord- 
ingly. 

The derivation of the QMC algorithm in Sec. lA II is 
valid for any dimension D of the lattice under consider- 
ation. Nevertheless, here we only report results for the 
sign problem in D = 1, the case considered in Sec. IIVI 
and make some remarks about the influence of the di- 
mensionality at the end. 

The average sign is defined as 



FIG. 12: Average sign (sign) [Eq. lIXM l of the fer mionic 
weight W{ as a function of electron-phonon coupling A in one 
dimension (a) for different inverse temperatures /3, and (b) 
for different values of the adiabatic ratio ljq. Here and in 
subsequent figures Ar = 0.05, lines are guides to the eye 
only, and errorbars are suppressed if smaller than the symbols 
shown. 



The inversion of the above relation is ill-conditioned and 
requires the use of the MEM.'^^ Fourier transformation 
leads to 



(sign) 



f/b 



(A20) 



The allowed imaginary times are = ^Ar, with nonneg- 
ative integers < ^ < i. Within the QMC approach, we 
have^ 

G\^{t,) = {G'^B,---B^)^^. (A21) 

Another interesting quantity is the one-electron density 
of states 

p(w - m) = - -Im G(w - /i) , (A22) 

TT 

where G{uj — /i) = N^^ G(fc, lo — ji). It may be ob- 
tained numerically via 



and subsequent use of the MEM. 



(A23) 



3. Sign problem 

Since a possible sign problem crucially affects QMC 
simulations, this section is devoted to a detailed inves- 



(kfl>t 



(A24) 



with the fermionic weight in the present case given by 
uif = det(l -I- Vt) [Eq. (|A5|) ]. We begin with the depen- 
dence of (sign) on the electron-phonon coupling strength. 
For simplicity, we show results for n = 0.5, while the ef- 
fect of band filling will be discussed later. The choice 
n = 0.5 is convenient since we know the chemical poten- 
tial, and we shall see below that the sign problem is most 
pronounced for a half-filled band. Moreover, all existing 
QMC results for the spinless Holstein model are for half 
filling, so that it is interesting to see how the sign problem 
affects simulations within the current approach. 

Figure I12r a) reveals that the average sign takes on 
a minimum near the point where — depending on the 
filling — the cross over to small polarons or a CDW state 
occurs. In the adiabatic regime cJq < 1, the cross over 
condition is A > 1 [Fig. I12r a')]. whereas for ZUq > 1 we 
have 5 > 1 [Fig. I12r b'l]. At weak and strong coupling, 
the average sign is close to 1, so that accurate simulations 
can be carried out. As shown in Fig. I12r b). we find that 
(sign) becomes very small for ZUq < 1 and intermediate A, 
while it increases noticeably in the nonadiabatic regime 
(Ug > 1. Owing to the absence of any autocorrelations in 
our approach, the present method therefore represents a 
significant improvement of existing algorithms, as the lat- 
ter face severe autocorrelations even for tUg > 1, thereby 
limiting the temperature range and cluster size. 

As illustrated in Fig. I13f a). the average sign depends 
strongly on the band filling n. While (sign) « 1 in the 
vicinity of n = or n = 1, a significant reduction is vis- 
ible near half filling n — 0.5. The minimum occurs at 
n — 0.5, and the results display particle-hole symmetry 
as expected. Here we have chosen j3t = 8, oJg = 0.4 and 
A = 1, for which the sign problem is most noticeable 
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FIG. 13: Average sign as a function of (a) band filling n, and 
(b) system size A'". 

according to Fig. ^1 Note that there is no sign prob- 
lem in simulations for, e.g., the Hubbard model at half 
filling, where the (fermionic) weights of the f and | elec- 
trons are equal and their product — entering the partition 
function — is positive.'*'' This can also be expected within 
the current approach if we consider the half-filled spin- 
full Holstein or Holstein- Hubbard model, and work along 
this line is in progress. 
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